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ABSTRACT 

In this paper we present a multiresolution-based method for period determination that 
is able to deal with unevenly sampled data. This method allows us to detect super- 
imposed periodic signals with lower signal-to- noise ratios than in classical methods. 
This multiresolution-based method is a generalization of the wavelet-based method 
for period det ection that is unable to deal with unevenly sampled data, presented by 
the authors in lOtazu et al.l l)2002j) . This new method is a useful tool for the analysis 
of real data and, in particular, for the analysis of astronomical data. 
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1 INTRODUCTION 

The detection of periodic signals in astronomical data 
■ has been usually addressed by classical Fourier-based or 
epoch folding methods. These methods have different prob- 
lems when dealing with non-sinusoidal periodic signals 
or with very low signal-to-noise ratios. When the anal- 
ysed data set contains several periodic signals, the behav- 
ior of classical period determination methods highly de- 
pends on intrinsic signal characteristics (see, for example, 
[ Lafler fc Kinmanl Il965t Ijurkevich]ll97lt IStellingwerJll978t 
| Lomblll976l: IScargld Il98l Il989l iRoberts. Lehar fc Dreherl 
Il987t iPress fc Rvbickvlll989l) . The reader is referred to the 
introduction of lotozuetalT (l2002t) (hereafter Paper I) for 
a general discussion. To avoid this problem, we presented 
there, as a preliminary step, a wavelet-based approach that 
only works with evenly time sampled data. 

In astronomy, however, this is not a usual situation, 
since data are mostly acquired on irregular intervals of time. 
In such a case there are two possibilities: resample the data 
into a new evenly sampled data set, or use a method able to 
deal with the original unevenly sampled data set. In the first 
case we are forced to modify the original data, which nec- 
essarily implies a loss of information. Moreover, this is not 
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always possible if the temporal gaps are larger than some 
of the periods present in the data. In order to avoid these 
problems, a technique capable to deal with unevenly sam- 
pled data is needed. 

The present paper is a natural extension of Paper I that 
allows to work on unevenly time sampled data. We show how 
the methodology of multiresolution decomposition (similar 
to the wavelet decomposition philosophy) is very well suited 
to this problem, since it is completely oriented towards de- 
composing functions into several frequential characteristics. 

As in Paper I, the main objective is to isolate every 
signal present in our data and to analyse them separately, 
avoiding their mutual influences. In Section [5] we outline 
some concepts in multiresolution analysis and their simi- 
larities with wavelet theory that are relevant to the stated 
problem. In Section[3]we propose an algorithm to detect each 
of the periodic signals present in a data set by combining 
multiresolution analysis decomposition with classical period 
determination methods. In Sections[I]and[S]we present some 
examples of synthetic data we used to test the algorithm and 
the results we obtained. We summarise our conclusions in 
Section H3 



2 MULTIRESOLUTION ANALYSIS 

Multiresolution decomposition introduces the con- 
cept of the presence of details between successive 
levels of resolution. Many wavelet decomposition 
algorithms are based on multiresolution analysis 



2 X. Otazu et al. 



schemes l|Chui| Il992l iDaubechies 




iMever 




199S ; 


Kaiser 1994; Vetterli & Kovaccvic 


1995; 


Mallat 


1999; 


Starck & Murtaebl 120021). and some astronomical amplica- 
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l200ll: IPolvgiannakis. Prcka-PaDadcma & Moussasl l200.tl 
In fact, all these methods use the same philosophy and 
the obtained results can be interpreted in the same way. 
However, wavelet theory presents some constraints on math- 
ematical functions. These constraints are not respected 
in all multiresolution decomposition schemes. Therefore, 
when we are using a certain multiresolution decomposition 
algorithm that fulfills the wavelet constraints, we are 
obtaining a wavelet decomposition. In contrast, when a 
given multiresolution decomposition algorithm violates 
these wavelet constraints, we are obtaining a result similar 
to, but which is not, a wavelet decomposition. 

As discussed in Paper I, we note that there are wavelet 
approaches that are based on approximations of a contin- 
uous wavelet transform and on th e subsequent study o f 
wavelet space coefficients (see, e.g., ISzatmarv et al.lll994t) . 
which are able to deal with unevenly sampled data sets. 
However, these algorithms present a non-direct inverse 
wavelet transform, in the sense that the search for period- 
icities is based on the fit between the wavelet base function 
profile and the signal one, and therefore on the values of 
the wavelet transform coefficients, which highly depend on 
the wavelet base used. Moreover, the period analysis has 
to be performed on the wavelet coefficients space but, since 
it is usually decimated, accurate period detection is a diffi- 
cult task. The multiresolution decomposition scheme we use 
in this work performs the decomposition on the temporal 
space, which allows to find accurate values for periodicities. 



tion frequency planes and p r is a residual signal (in fact 
n — r, but we explicitly substitute n by r to clearly express 
the concept of residual). In our case, we are using a dyadic 
decomposition scheme. This means that the standard devi- 
ation of the Gaussian function associated with the Fi filter 
is o~i = 2oi-\. Thus, similarly to the wavelet approach, the 
original signal po has double resolution compared to pi , and 
so on. All these pi (i = 0, ■ ■ ■ , n) signals have the same num- 
ber of data points as the original signal p. 

Since ct;+i depends on en, a value for ao has to be care- 
fully chosen for every data set. It has to be fixed considering 
a likely minimum value for the time-duration of the features 
and the characteristic time sampling of the data set, in order 
to include a significant number of points on which to per- 
form the convolutions. Too small values do not accurately 
describe feature profiles, and suffer from poor or noisy data. 
In contrast, too large values reduce the noise effect by inte- 
grating a lot of data points, but may ignore interesting high 
frequency features. Hence, a first analysis of data has to be 
performed in order to obtain a useful initial ao value. 

We have used the same notation as in the wavelet 
decomposition described in Paper I because, as explained 
above, the idea of these multiresolution planes is similar to 
the wavelet ones. 

We note that this particular decomposition scheme that 
uses a Gaussia n kernel can also be interpreted as a scale- 
space filterin g teaubaud et alJll986t ISoorring et"al"lll997l . 
IWitkinl lT994) or as a pa rticular case of more general im- 
age d iffusion approaches (lLindeberdll990| . iPerona fc Malikl 
1990). Smoothed data sets pi can be interpreted as diffused 
scale-space images, and the difference between them as the 
details at different scales. 



2.1 A multiresolution analysis algorithm 

In order to obtain a multiresolution decomposition for sig- 
nals, an algorithm to decompose the signal into frequency 
planes can be defined as follows. Given a signal p we con- 
struct the sequence of approximations: 

pl=Fl(p), P2=F' 2 (pi), P3 = F 3 (p2), ■ ■ ■ , (1) 

performing successive convolutions with Gaussian filters Fi. 

It is important to note the difference between this se- 
quence of convolutions and the one used in Paper I. In the 
latter we were dealing with a discrete convolution mask, and 
hence forced to work with evenly spaced data. In contrast, 
the continuous nature of the convolution functions used in 
the present paper, allows us to work with unevenly sampled 
data. 

Similarly to the wavelet planes, the multiresolution 
frequency planes are computed as the differences between 
two consecutive approximations and pi. Letting Wi — 
Pi-i — pi (i = 1, • • • , n), in which po = p, we can write the 
reconstruction formula: 
n 

P = ^ Wi + p r . (2) 

; = 1 

In this representation, the signals pi (i = 0, • ■ • , n) are 
versions of the original signal p at increasing scales (decreas- 
ing resolution levels), Wi (i =!,■■■ ,ri) are the multiresolu- 



3 PERIOD DETECTION ALGORITHM 

We propose to apply this multiresolution analysis algorithm 
to solve our initially stated problem: to isolate each of the 
periodic signals contained in a set of unevenly sampled data 
and study them separately. 

In order to do so, we proceed as in Paper I, and the 
period detection algorithm we propose is as follows: 

(i) choose values for ao and n; decompose the original 
signal p into its multiresolution frequency planes Wi (i — 
0, •••,«). 

(ii) detect periods in each of the n obtained frequency 
planes Wi. 

P hase Dispers ion Minimization ( pdm) iStellingwerj 
Il978l) and CLEAN iRoberts et al.lll987fl methods are used 
to detect periods in the original data and in every mul- 
tiresolution frequency plane. In Paper I we described sev- 
eral undesirable effects of PDM and CLEAN methods on data 
with superimposed signals, as well as the advantages of using 
multiresolution-based methods over the classical ones. 

Hereafter, and for notational convenience, the 
multiresolution-based PDM and CLEAN methods will 
be called MRPDM and MRCLEAN, respectively. 
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4 SIMULATED DATA 

In order to check the benefit of applying MRPDM versus 
PDM, or MRCLEAN versus CLEAN, we proceeded similarly as 
in Paper I by generating and analysing several sets of sim- 
ulated data containing two superimposed periodic signals. 
Each data set is composed of a high-amplitude primary si- 
nusoidal function and a secondary low-amplitude Gaussian 
one. Finally, we added a white-Gaussian noise to this com- 
bination of signals. We increased the value of the noise stan- 
dard deviation, a, up to the value where detection of periodic 
signals became statistically insignificant in both the classical 
and the multiresolution-based methods. 

The primary signal is intended to simulate variable 
sources with a pure sinusoidal intensity profile (like preces- 
sion of accretion discs), and the secondary, burst-like events 
(like pulses or eclipses) superimposed to it. We note that in 
Paper I we also studied the superposition of two sinusoidal 
functions, the so called sine+sine case. However, here we 
have directly focused on the more realistic situation of the 
so called sine+ Gaussian case, where the difference between 
the classical and the multiresolution-based methods is more 
critical (see Paper I). 

The characteristics of the signals are: 

(i) Each signal is generated as an unevenly spaced data 
set. Its time sampling has been taken from observations of 
the X-ray binary LM C X-l by the All Sk y Monitor onboard 
the RXTE satellite llLevine et alJ 1996). We have used an 
observation period which has 8270 measurements during 679 
days. 

(ii) The high-amplitude sinusoidal function has an ampli- 
tude equal to 1 and a period of 108.5 days. 

(iii) The amplitudes of the low-amplitude periodic Gaus- 
sian function are 0.1 and 0.5. 

(iv) The periods used for the secondary function are 13.13 
and 23.11 days. 

(v) We have used two values for the full width at half- 
maximum (FWHM) of the Gaussian signal, corresponding 
to 2 and 6 days, respectively. 

In the first four columns of Table Q we present the pa- 
rameters used to generate each simulated data set. 



5 RESULTS 

We recall that the simultaneous use of two independent 
methods, such as CLEAN and PDM, is usually applied to 
discriminate false period detections from the true ones. A 
similar procedure can be used with each one of the mul- 
tiresolution planes in the MRPDM and MRCLEAN methods. 
Therefore, when comparing the behavior of these different 
methods, we have to compare the usual PDM-CLEAN method 
combination for period estimation prior to the new MRPDM- 
MRCLEAN combination. 

Taking into account the characteristics of the simulated 
data we have chosen a <to = 1 day. As in Paper I, the pri- 
mary period (108.5 days) is always detected by all methods, 
and does not appear in Tabled In this table, the four last 
columns show the detected low-amplitude periods for each 
data set using PDM, clean, mrpdm and mrclean meth- 
ods, respectively. A dash is shown when a period is not de- 
tected, and a question mark when the detection is difficult 



or doubtful. When a period is found in the multiresolution- 
based methods, we also show in parentheses the mutiresolu- 
tion planes where it is detected. 

We must note that the use of two different FWHM for 
the Gaussian, combined with two different periods (13.13 
and 23.11 days), gives 4 different profiles. Hence, the phase 
duration of the burst-like event ranges from very low to rel- 
atively high values in the following order: FWHM=2 and 
period=23.11, FWHM=2 and period=13.13, FWHM=6 and 
period=23.11, and finally FWHM=6 and period=13.13. 

In view of the results displayed in Tabled we can make 
the following comments: 

(i) In all methods, with high noise-to-signal ratios the de- 
tected periods are slightly different from the simulated ones. 

(ii) There is a better performance of CLEAN over PDM. 
We must note that when the FWHM is only 2 days, PDM 
never detects the secondary period. Only with FWHM=6 
days and a relatively high amplitude (0.5), can PDM detect 
the low-amplitude periodic signals. 

(iii) As the noise increases, the detection starts to fail 
in the lower multiresolution planes (higher frequencies), 
and only the higher ones (lower frequencies) are noise-free 
enough to allow period detection. 

(iv) MRPDM and MRCLEAN perform better or similar 
than PDM and CLEAN methods (see exception below). When 
CLEAN marginally detects the secondary period, MRPDM 
and most of times MRCLEAN have no problems to detect 
it, and they work properly even with higher noise. In the 
MRPDM case, the results are always better than with PDM. 

(v) In all cases with FWHM=2, the MRPDM performance 
is much better than MRCLEAN, because the signal is clearly 
non-sinusoidal. In the FWHM=6 cases, MRPDM is only 
slightly better than MRCLEAN, since the signals are closer 
to a sinusoidal profile. 

(vi) For a given amplitude of the Gaussian signal, the 
maximum noise-to-signal ratio achieved with MRPDM and 
MRCLEAN increases with the phase duration of the FWHM. 

The only exception to |(iv)| is for the most extreme of 
the simulated cases, i.e., the one with the lower amplitude, 
lower FWHM value and longer period. However, we note 
that the period detected by CLEAN is slightly different than 
the simulated one. 

All these results are very similar to those shown in Ta- 
ble 2 of Paper I. Nevertheless the maximum noise-to-signal 
ratios achieved in the present cases are around 2.5 times 
higher. This can be explained because, although the time 
span of the data sets used here is around 1.5 times smaller, 
the number of points per unit time is around 12 times higher 
than in Paper I. 

Finally, and for illustrative purposes, we show in Fig. Q 
the simulated data set generated with the following: 13.13- 
day period, FWHM=6.0 days and amplitude=0.1 with 
Onoise = 0.6. The outputs of PDM and CLEAN, are also 
shown. None of these methods is able to detect the 13.13-day 
period. We show in Fig. [5] the outputs from MRPDM. The 
simulated period is detected, with its corresponding sub- 
harmonics, in the multiresolution planes W2 and W3. The 
MRCLEAN outputs, shown in Fig. [3J reveal a marginal de- 
tection of the 13.13-day period in multiresolution plane W3. 
We note that we would not consider this detection as a true 
one when taken alone. However, since the same period is 
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Table 1. Periods, in days, detected in the data sets. The first four columns are the parameters used 
for the simulated data. A dash is shown when a period is not detected, and a question mark when the 
detection is difficult or doubtful. When a period is indicated in the multiresolution-based methods, we 
also show in parentheses the multircsolution planes where it is detected (using <r = 1 day). 
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SIMULATED DATA 
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Figure 3. MRCLEAN output for each wavelet plane of the data 
at the top of Fig. 



Figure 1. Top: simulated data of a Gaussian with 13.13 days 
period, a FWHM of 6.0 days, amplitude=0.1 and rj no i se = 0.6, 
superimposed on the sinusoid with 108.5 days period and ampli- 
tude^. 0. Bottom: PDM and CLEAN outputs of this data set. 
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Figure 2. MRPDM output for each wavelet plane of the data at 
the top of Fig. [I] 



clearly detected in two multiresolution planes of MRPDM, 
we can establish the existence of this period in the analysed 
data set. 



6 CONCLUSIONS 

In this paper we have presented a multiresolution-based 
method for period determination able to deal with unevenly 
sampled data. This constitutes a significant improvement 
with respect to the wavelet-based method presented in Pa- 
per I, which is unable to deal with unevenly sampled data. 
The overall performance of the present method is similar to 
the wavelet-based one, in the sense that it allows us to detect 
superimposed periodic signals with lower signal-to-noise ra- 
tios than in classical methods. We stress that one advantage 
of the present method over classical methods is the simulta- 
neous detection of a period in more than one multiresolution 
plane, allowing to improve the confidence of a given detec- 
tion. Moreover, since here we are not forced to lose or modify 



the information when averaging or interpolating the original 
data, we can reach higher noise-to-signal ratios than in the 
wavelet-based method described in Paper I. 

We note that the multiresolution decomposition scheme 
that we have used can be interpreted as a particular case 
of scale-space filtering. In order to improve isolation of pe- 
riodic features, more general approaches could be used to 
perform this decomposition . In this context, anisot ropic dif- 
fusion schemes proposed bv lPerona fc Malikl il990f) could be 
useful if properly tuned. 
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